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I.  Introduction 

The  use  of  second  order  or  higher  interpolation  schemes  for 
convection  terms  causes  the  appearance  of  unwanted  extremes  near 
points  where  the  transported  function  changes  its  slope  rapidly.  The 
donor  cell  method  uses  first  order  interpolation  and  therefore  cannot 
generate  such  extrema,  but  it  has  a diffusion  which  is  unacceptable 
in  most  cases.  The  partial  donor  cell  method  (PDM)  presented  here  is 
nonlinear  combination  of  a higher  order  scheme  and  the  donor  cell 
method.  It  adds  just  enough  diffusion  to  prevent  the  occurrence  of 
such  extrema.  As  compared  with  other  hybrid  methods  it  is  less 
diffusive  because  it  restricts  the  additional  diffusion  to  points 
where  it  is  absolutely  necessary.  Flux-Corrected  Transport  (FCT) , on 
the  other  hand,  removes  the  extrema  after  they  have  occured.  The 
results  of  this  method  and  FCT  are  very  similar,  as  demonstrated  in 
the  test  runs.  The  PDM  leaves  the  solution  undisturbed  for  zero 
velocities.  It  therefore  cannot  be  used  to  create  artificial 
viscosity  in  stationary  or  nearly  stationary  shock  fronts.  The 

method  can  be  readily  extended  to  multidimensions. 

II.  Theory 

Let  f be  the  function  to  be  transported  across  a one-dimensional 
grid;  we  denote  by  f^  the  value  of  the  function  at  x ^ . Furtnermore, 
let  d be  the  difference  operator, 

‘*^j+l/2  “ ^j+1  “ ^j- 
Note:  Manuscript  submitted  January  30,  1978. 
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Two  higher-order  schemes  are  considered  here. 


a)  A simple  scheme. 

Let  all  quantities  be  defined  at  integral  points  at  space 
and  time,  with  the  exception  of  the  velocity  u,  defined  at  half 
integral  points  in  space  and  time.  Then,  the  new  value  of  f^  is 


given  by 


dt 


f j - f j - 2 ^“j+i/2  ^^j+1  ■ ^j-1/2  ^j-l^^^‘^^j+1/2 

The  scheme  is  unstable,  with  an  amplification  factor 

2 


A = 1 + 1/2 


The  addition  of  a second  order  term,  corresponding  to  a second  order 
interpolation,  will  make  this  scheme  stable, 
b)  Lax-Wendroff  Scheme. 

In  this  scheme,  one  defines  intermediate  values  as 

dt 


^j+1/2  ■ 2 ^‘"j+l  ^j+1  “ 

'V 

then,  with  defined  similarly 

A 

^‘^j+1/2  ^j+1/2  ■ '^j+1/2  ^j-l/2^/‘^*j 
The  donor  cell  method  can  be  written  in  the  form 

'j  ■ - "‘8»(Uj))  f u -£j„j)/dx.^^^2 

(l+sl8„(Uj))  - £j.i 

or 


(3) 


(4) 


(5) 


2 


[^^i+1  “j+1  ■ - ^j-l“j-l>/‘^*j-l/2] 

dt  r 

+ ^ sign(Uj)^(fj^^  - fjU^)/dx.^^/2  - - ^n-l>/‘^*j-l/2j  ’ 


As  is  well-known,  the  difference  between  the  donor  cell  method  and 
* 

the  simple  scheme  (a)  defined  above  is  a diffusion  term  (the  second 
term  in  equation  5a)  first  order  in  time  and  second  order  in  space. 
(The  difference  in  the  case  of  the  Lax-Wendroff  scheme  also  involves 
a higher  order  term  in  space  and  time.) 

In  order  to  derive  the  partial  donor  cell  method,  let  us  define 
for  the  scheme  (a) 

^j+1/2  ° 1 ''j+1/2  I ^ ‘^*j+l/2,  (6) 


and  for  scheme  (b) 

^j+1/2  “ I I '^‘^^j+1/2  “ ‘^’^l“j+l/2l^‘^*j+l/2) 


Furthermore,  let  a 


j+1/2 


be  an  array  of  numbers  such  that 


° - °j+l/2 


< 1; 


here  o represents  the  fraction  of  diffusion  to  be  added.  If  a = 1, 
the  complete  donor  cell  method  would  result.  It  should  be  emphasized, 
however,  that  this  holds  for  variable  velocities  only  up  to  terms 
second  order  in  time  and  third  order  in  space. 

With  these  definitions,  we  add  the  following  diffusion  terms  in 
a conservative  manner  to  either  the  simple  scheme  or  the  Lax-Wendroff 
scheme: 


° ^j+1/2  °j+l/2  ‘*^j+l/2  " ^j-1/2  °j-l/2  ‘^^j-1/2 
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In  order  to  simplify  the  following  discussion  let  us  assume  for  a 
moment  constant  velocity  and  u > 0,  and  discuss  only  the  simple 
scheme.  Then  with 

e = u dt/dx 

we  can  rewrite  the  transport  algorithm  including  diffusion  as 

~ I ^^j+1  " ^j-1^  2 ^°j+l/2  ‘^^j+1/2  ■ °j-l/2  ‘^^j-1/2^’ 


or  setting 


‘'3+1/2  °j+l/2  ‘^^j+l/Z, 


- 2 ^^j+l  ■ 2 ^‘'j+1/2  ~ ‘'j-1/2^* 


(9) 


(10) 


J J 

The  main  question  is  how  to  determine  p.  As  the  convection  is  nothing 
but  an  interpolation  between  and  (for  u > 0) , the  new  value 

of  the  function  must  lie  between  f^  and  This  implies  if  f^  is 

not  an  extrema  in  the  absence  of  diffusion  that 

1 1 Vi- Vi  I I'j  - Vil- 

If  this  is  not  the  case,  some  diffusion  has  to  be  added  to  make  the 
updated  function  value  lie  between  these  limits.  The  diffusion 
depends  obviously  on  the  direction  in  which  the  fluid  moves  and 
therefore  should  incorporate  this  information.  One  way  to  take  it 
into  account  is  the  following.  Let  us  define 

Sj  = A + I I sign(dfj^^/2>  + sign(dfj_j^^2>  I’  (11) 

then  we  write  general  in  the  form 
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^j+1/2  = (df  >‘^^j+l/2l 

1/2  (1  - sign  1 ‘*^j+3/2l 

- 1/2  Sj  (1  + ,ig„  I ■‘*1-1/2  I] 


(12) 


This  prescription  means  we  subtract  df,  the  difference  dfj^^y2  taken 
in  the  direction  of  the  flow,  the  amount  A whether  or  not  the  func- 
tion has  an  extremum  at  this  point,  and  the  amount  B only  if  there  is 
no  extremum.  But,  if  the  amount  exceeds  the  difference 
diffusion  will  result.  This  is  the  point  where  (as  in  FCT)  the  non- 
linearity comes  in.  The  values  of  A and  B remain  to  be  determined. 

In  order  to  relate  this  diffusion  to  the  hybrid  schemes 
2 3 

described  by  Harten  and  Van  Leer  we  define  a quantity  0 by 


l‘^^j+l/2l  ■ l^^j-l/2l 


- 


|df, 


(13) 


j+l/2l  l‘^^j-l/2l 


(compare  Eq.  4.12a,  4.12b  of  ref.  2).  Marten's  method  then  can  be 
written  in  our  notation  as 
H 


^.j+1/2  ®j+l^  ‘^^j+1/2 


(14) 


This  means  the  solution  has  an  added  diffusion  proportional  to  this 
smoothness  parameter  6^ . In  order  to  make  the  connection  with  our 
methods  more  obvious,  let  us  define 

0j  “ l‘*^j+l/2l  ■ l‘^^j-l/2l^^l‘*^j+l/2l’ 

0-  - max  (0,  |dfj_^/2l  " 1^^  j-n/zl ) / j.1/2 1 • ^^b) 
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With  these  definitions  we  can  rewrite  Eq.  (12)  as 

^j+1/2  ° ‘^^j+1/2  ^“j+1/2^^  ®j+l 

+ (1  + sign 

The  choice  of  A = 1,  B = 0,  therefore,  closely  resembles  the  hybrid 
scheme  of  Harten  and  Van  Leer. 

The  important  difference  is  that  in  the  partial  donor  cell 
method  only  as  much  diffusion  is  added  as  is  needed  to  keep  the 
function  f monotone.  This  is  achieved  in  two  ways;  first,  by  consid- 
ering only  the  differences  in  the  direction  from  which  the  fluid 
comes;  second,  by  taking  the  maximum  in  Eqs.  (15).  The  possibility 
of  adjusting  the  parameters  A and  B is  a further  advantage. 

We  now  discuss  the  manner  in  which  A and  B are  determined.  The 
simplest  argument,  which  leads  to  A = 0 and  B = 1,  is  the  following. 
If 

1/2  I I i I tj  - fj.i  I 

then  no  extrema  can  occur.  That  is,  if  f^  has  no  extrema,  this  is 
equivalent  to 

1 f . ,1  - f . 1 < I f . -f . J 

Therefore,  no  additional  diffusion  is  needed  implying  A + B = 1. 

The  choice  of  linear  interpolation  at  all  extrema  leads  to  A = 0.  A 
more  sophisticated  argument  proceeds  as  follows.  It  is  presented 
here  for  scheme  a:  extreme  and  as  long  as 


2 I ^j+1  ■ ^j-1  I ^ ' ^j  ■ ^j-1'’ 
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[f  one  chooses  B such  that  only  additional  diffusion  is  added  if 
this  condition  is  not  satisfied,  namely  (similar  for  b)) 

a)  B = - - 2,  b)  B = - 

e e 

then  at  points  where  the  subtracted  difference  is  smaller,  f.  = f.  , 

J .1-1 

regardless  of  the  value  of  e.  This  leads  to  great  phase  errors 
locally.  A good  choice  is  to  take  the  values  found  for  e = 1/2.  The 
results  for  the  two  schemes  are 

a)  A = 1 b)  A = 1 

B = 2 B = 4 

(One  should  remark  that  usually  in  hydrodynamics  e at  its  maximum  is 
about  0.3;  therefore,  the  phase  errors  should  not  be  too  big.)  This 
choice  assures  the  correct  range  of  f^.  The  test  runs  confirm  also 
that  it  is  less  diffusive  than  the  choice  A = 0,  B = 1.  The  results 
are  for  all  practical  purposes  identical  with  those  of  FCT  for  the 
transport  of  square  waves.  Test  runs  with  rapid  varying  velocity 
have  shown  that  the  best  results  were  achieved  with  A = 0,  B = 1.  A 
description  of  this  algorithm  is  given  in  Appendix  B. 

Test  runs; 

Test  runs  have  been  done  for  three  cases: 
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1)  Squarewave  10  meshpolnts  across. 

2)  Squarewave  two  meshpolnts  across. 

3)  Cosine  wave  with  10  points  for  the  total  period. 


For  the  choice  e = .2,  the  results  are  presented  at  10  and  100  steps, 
for  which  the  solution  is  propagated  2 and  20  mes..i'Oints,  respec- 
tively. The  results  are  given  in  the  tables  below  (tables  were 
chosen  because  of  the  small  differences).  One  sees  clearly  that  for 
all  the  test  A = 0,  B =■  1,  is  most  diffusive,  but  gives  basically 
the  same  results.  The  choice  of  A = 1,  B = 2 in  the  10  point  sqiare 
wave  yields  the  same  results  as  does  FCT  (as  described  in  the 
appendix) . The  results  of  the  other  tests  seem  to  show  that  the  PDM 
holds  the  extrema  a little  bit  better  than  FCT.  This  is  clearly 
demonstrated  in  the  cosine  wave  case  where  the  maximum  after  100  steps 
is  1.729  compared  to  1.518  for  FCT.  Also,  the  typical  three-point 
plateau  for  FCT  appears,  whereas  PDM  keeps  the  maximum  at  the  right 
place.  Each  of  the  schemes  has  a tendency  to  form  a plateau  with 
more  extension  in  the  direction  of  the  flow  for  a maximum  and  vice 
versa  for  a minima.  The  choice  of  A = 0,  B = 1 for  the  simple 
scheme  seems  to  have  this  tendency  in  far  lesser  degree.  After 
transport  over  20  meshpolnts  the  cosine  wave  is  practically  symmet- 
rical, but  with  about  half  the  amplitude  found  for  the  Lax-Wendroff 
scheme. 

III.  Conclusions 

A simple  transport  scheme  has  been  demonstrated.  It  prevents 
the  occurrence  of  new  extrema  (and  therefore  assures  monatonicity) . I 

I 

i 
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The  additional  donor  cell  contributions  needed  can  be  chosen 
according  to  the  two  parameters  A and  B.  The  most  diffusive  choice 
A = 0,  B = 1,  will  result  in  a scheme  slightly  more  diffuse  than  FCT, 
but  it  will  guarantee  the  nonoccurrence  of  new  extrema  caused  by  the 
higher  order  scheme.  Application  of  PDM  in  one  dimension  requires 
the  same  number  of  operations  as  FCT  (for  the  ASC  at  NRL) . It  can 
readily  be  applied  to  nonsplit  multidimensional  problems. 
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Appendix  A 


In  this  appendix  a short  description  of  the  FCT  algorithm 
actually  used  in  the  test  runs  is  given.  Theory  and  more  details 
are  given  in  ref.  (1).  Let  H be  the  transport  operator  for  one 
timestep  (either  the  simple  scheme  (a)  or  the  Lax-Wendroff 
scheme) , 
then 

= Hf . (Al) 

if  is  the  transported  function.  Now  a diffused  function  f^^  is 
computed  using  the  old  values 


f"  -f  1/8 


j+1/2 


(A2) 

, it  is  useful  to 


in  order  to  obtain  the  antidif fusive  flux  df 
define  a factor  S by 

^j+1/2  ° ^‘^^j+1/2^ 

Cisign  (df“3/2)  + sign  (df°!J3/2>  + (df (A3) 

which  is  0 if  the  sign  of  the  three  differences  are  different*  With 
the  help  of  this  factor,  the  antidif fusive  flux  is  given  by 

‘^^j+l/2  ° ^j+1/2  l‘^^j-l/2l 


Then  final  resulting  function  is  then  given  by 


/jfA  jfA  \ 

^‘^^f+1/2  ■ '*^j-l/2^ 


(A5) 
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Appendix  B 


Here  a short  description  for  the  PDM  scheme  for  variable 
velocities  and  grids  are  given.  The  choice  of  A = 0,  B =1  is 
incorporated.  Results  obtained  in  a multidimensional  hydrocode 
applying  PDM  have  been  optimal  with  this  choice. 

Let  H be  the  hydrooperator  for  one  timestep 

f^  = Hf. 
j J 


(Bl) 


then  fj  is  the  updated  function. 

The  following  formula  can  be  applied  in  different  directions  in 
a onestep  hydrocode  by  taking  the  appropriate  differences.  Be 
Sj  = 1/4  I sign  (dfj^j^/2^  ^‘^^j-1/2^  I 

then  Sj  vanishes  if  the  signs  are  different.  It  is  1/2  for  the  same 
signs.  The  diffusive  flux  difj^j^^2)  then  given  by 

‘^^^j+1/2  “ dt  ^‘^^j+l/2>  l‘^fj+l/2l 


- (1  - sign  |dfj^^^2l 

- (X  + sign  l‘^^j-l/2l 

Finally  the  new  values  f are  given  by 


fj  = fj  + (difj^j^^2  ‘^^^j-l/2^^‘**j+l/2 


(B3) 


(B4) 
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Analytic 

Values 


t = 2.0 


t = 20.0 


.012 

.032 

0.0 

0.000 

.000 

.000 

.002 

.077 

.000 

0.0 

.006 

.039 

.004 

.040 

.169 

.004 

1 

0.0 

.247 

.264 

.246 

.409 

.339 

.436 

1.0 

.789 

.745 

.795 

.733 

.620 

.739 

! 

1.0 

.963 

.758 

.964 

.886 

.847 

.889 

l.Q 

.995 

.996 

.996 

.955 

.942 

.956 

1.0 

1.000 

1.000 

1.000 

.984 

.974 

.984 

1.0 

1.000 

1.000 

1.000 

.995 

.980 

.995 

1.0 

1.000 

1.000 

1.000 

.998 

.978 

.998 

1.0 

1.000 

1.000 

1.000 

.998 

.964 

.998 

1.0 

1.000 

1.000 

1.000 

.997 

.924 

.998 

1.0 

.994 

.962 

1.000 

.959 

.834 

.996 

1.0 

.753 

.736 

.754 

.592 

.663 

.565 

0.0 

.211 

.255 

.206 

.268 

.382 

.2bl 

0.0 

.037 

.042 

.036 

.114 

.159 

.111 

0.0 

.005 

.005 

.005 

.045 

.062 

.044 

.016 

.023 

.016 

.006 

.008 

.006 

.002 

.002 

.002 

TABLE  1. 

10  Point  Square  Wave 
Simple  Method 
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TABLE  2 

2-Point  Square  Wave 
Simple  Method 


I- 

C 
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.000 

.069 

.132 

.148 

.274 

.547 

.482 

.191 

.125 

.237 

.148 

.326 

.639 

.482 

.691 

.666 

.612 

.604 

.940 

.816 

.849 

1.309 

1.460 

1.345 

1.485 

1.504 

1.180 

1.363 

1.809 

1.864 

1.814 

1.810 

1.700 

1.381 

1.503 

2.000 

1.931 

1.868 

1.852 

1.726 

1.403 

1.518 

1.809 

1.875 

1.763 

1.852 

1.674 

1.361 

1.518 

1.309 

1.334 

1.389 

1.396 

1.060 

1.184 

1.151 

.691 

.541 

.655 

.515 

.496 

.820 

.oj7 

.191 

.136 

.186 

.191 

.300 

.619 

.497 

.000 

.069 

.132 

.148 

.274 

.597 

.482 

TABLE  3 


1 


Simple  Method 
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TABLE  4. 


10-Point  Square  Wave 


LW  Method 


t = 

Analytic 

Values 

2.0 

t = 20.0 

A 

1.0 

0.0 

FCT 

A 

1.0 

0.0 

FCT 

B 

4.0 

1.0 

B 

4.0 

1.0 

0.0 

.000 

.000 

.000 

.000 

.028 

.000 

0.0 

.000 

.000 

.000 

.009 

.061 

.010 

0.0 

.000 

.000 

.000 

.153 

.122 

.138 

0.0 

.003 

.047 

.006 

.298 

.215 

.248 

0.0 

.390 

.335 

.432 

.368 

.282 

.303 

1.0 

.762 

.659 

.576 

.378 

.300 

.314 

1.0 

.579 

.604 

.576 

.331 

.296 

.314 

0.0 

.210 

.285 

.332 

.214 

.268 

.300 

0.0 

.048 

.060 

.069 

.125 

.203 

.186 

0.0 

.008 

.009 

.010 

.067 

.114 

.000 

0.0 

.001 

.001 

.001 

.033 

.054 

.050 

0.0 

.000 

.000 

.000 

.015 

.024 

.023 

TABLE  5. 

2-Point  Square  Wave 
LW  Method 
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t 


2.0 


t - 20.0 


Analytic 

Values 


A 

1.0 

0.0 

B 

4.0 

1.0 

A 

1.0 

0.0 

B 

4.0 

1.0 

.000 

.085 

.157 

.188 

.409 

.712 

.665 

.191 

.175 

.271 

.188 

.579 

.753 

.665 

.691 

.743 

.655 

.658 

1.067 

.894 

.939 

1.309 

1.391 

1.323 

1.429 

1.410 

1.135 

1.228 

1.809 

1.811 

1.772 

1.762 

1.571 

1.270 

1.325 

2.000 

1.915 

1.843 

1.812 

1.591 

1.288 

1.335 

1.809 

1.825 

1.728 

1.812 

1.421 

1.247 

1.335 

1.309 

1.257 

1.345 

1.342 

.934 

1.105 

1.061 

.691 

.603 

.677 

.571 

.590 

.865 

.772 

.191 

.189 

.229 

.238 

.429 

.730 

.674 

.000 

.085 

.157 

.188 

.409 

.712 

.665 

TABLE  6. 

Lax-Wendroff  Method 
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